# delimit ;
clear ;
cd "replication" ;
set more off ; 

* ***************************************************************************** ;
* multiple inference adjustment
* ***************************************************************************** ;

* first, we need to collect all the estimates from the main tables ; 
* make dataset to keep the estimates  ;

clear ; 
set obs 1000 ;

gen table = "" ;
gen outcome = "" ;
gen coefficient = "" ;
gen estimate = . ; 
gen pval = . ;

local i = 1 ;
save temp, replace ; 

* ***************************************************************************** ;
* collect results from midline knowledge 
* ***************************************************************************** ;

cd ".." ;
do "./replication/do/table03-midline-knowledge-aggregate.do" ;

* regressions: agggregate scores ; 
foreach var of varlist score_tot score_cleanliness score_midwife score_condom { ; 
	areg `var' healthonly healthandpay , a(strata) robust ;
	
	preserve ;
	use temp, clear ;
	
	foreach coeff in healthonly healthandpay { ;
		replace table = "midline knowledge" in `i' ;
		replace outcome = "`var'" in `i' ;
		replace coefficient = "_b[`coeff']" in `i' ; 
		replace estimate = _b[`coeff'] in `i' ;
		replace pval = 2*ttail(e(df_r),abs(_b[`coeff']/_se[`coeff'])) in `i' ;
		local ++ i ;
	} ;
	save , replace ;
	restore ;
} ;

* ***************************************************************************** ;
* collect results from endline knowledge 
* ***************************************************************************** ;

cd ".." ;
do "./replication/do/table04-endline-knowledge-aggregate.do" ;

foreach var in module_all module1_know module23_know module4_know module5_know { ;
	
	areg `var' healthonly healthandpay, a(strata) robust ;
	qui est sto `var'; 
		
	preserve ;
	use temp, clear ;
	
	foreach coeff in healthonly healthandpay { ;
		replace table = "endline knowledge" in `i' ;
		replace outcome = "`var'" in `i' ;
		replace coefficient = "_b[`coeff']" in `i' ; 
		replace estimate = _b[`coeff'] in `i' ;
		replace pval = 2*ttail(e(df_r),abs(_b[`coeff']/_se[`coeff'])) in `i' ;
		local ++ i ;
	} ;
	save , replace ;
	restore ;

} ;

* ***************************************************************************** ;
* collect results from endline behavior
* ***************************************************************************** ;

cd ".." ;
do "./replication/do/table05-endline-behavior.do" ;

foreach var in behav_all $module1_behav $module23_behav $module4_behav $module5_behav  { ;
	
	areg `var' healthonly healthandpay, a(strata) robust ;
		
	preserve ;
	use temp, clear ;
	
	foreach coeff in healthonly healthandpay { ;
		replace table = "endline behavior" in `i' ;
		replace outcome = "`var'" in `i' ;
		replace coefficient = "_b[`coeff']" in `i' ; 
		replace estimate = _b[`coeff'] in `i' ;
		replace pval = 2*ttail(e(df_r),abs(_b[`coeff']/_se[`coeff'])) in `i' ;
		local ++ i ;
	} ;
	save , replace ;
	restore ;
	
} ;

* ***************************************************************************** ;
* FDR adjustment of p-values
* ***************************************************************************** ;

use temp, clear ;

* remove empty rows ;
drop if missing(table) & missing(outcome) & missing(coefficient) & missing(estimate) & missing(pval) ;

* The following code is from Anderson (JASA, 2008) ; 
* Collect the total number of p-values tested ;
quietly sum pval ;
local totalpvals = r(N) ;

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank ;
quietly gen int original_sorting_order = _n ;
quietly sort pval ;
quietly gen int rank = _n if pval~=. ;

* Set the initial counter to 1 ;
local qval = 1 ;

* Generate the variable that will contain the BH (1995) q-values ;
gen bh95_qval = 1 if pval~=. ;

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, 
*then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses 
* are rejected at q = 0.998, etc. The loop ends by checking which hypotheses
* are rejected at q = 0.001. ;

while `qval' > 0 { ;

	* Generate value qr/M ;
	quietly gen fdr_temp = `qval'*rank/`totalpvals' ;
	
	* Generate binary variable checking condition p(r) <= qr/M ;
	quietly gen reject_temp = (fdr_temp>=pval) if fdr_temp~=. ;
	
	* Generate variable containing p-value ranks for all p-values that meet above condition ;
	quietly gen reject_rank = reject_temp*rank ;
	
	* Record the rank of the largest p-value that meets above condition ;
	quietly egen total_rejected = max(reject_rank) ;
	
	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition ;
	replace bh95_qval = `qval' if rank <= total_rejected & rank~=. ;
	
	* Reduce q by 0.001 and repeat loop ;
	quietly drop fdr_temp reject_temp reject_rank total_rejected ;
	local qval = `qval' - .001 ;
	
} ;
	
quietly sort original_sorting_order ;

* ***************************************************************************** ;
* make table for output, doing this manually 
* ***************************************************************************** ;

* make the columns ;
gen coeff = "" ; 
replace coeff = "\textit{Panel A. Effects on Short-Term Health Knowledge}" in 1 ;
replace coeff = "HEE, \(p\)-value" in 2 ; 
replace coeff = "HEE, \(q\)-value" in 3 ; 
replace coeff = "\\HEEC, \(p\)-value" in 4; 
replace coeff = "HEEC, \(q\)-value" in 5; 

replace coeff = "\\\textit{Panel B. Effects on Longer-Term Health Knowledge}" in 8 ;
replace coeff = "HEE, \(p\)-value" in 9 ; 
replace coeff = "HEE, \(q\)-value" in 10 ; 
replace coeff = "\\HEEC, \(p\)-value" in 11; 
replace coeff = "HEEC, \(q\)-value" in 12; 

replace coeff = "\\\textit{Panel C. Effects on Health Behavior}" in 13 ;
replace coeff = "HEE, \(p\)-value" in 14 ; 
replace coeff = "HEE, \(q\)-value" in 15 ; 
replace coeff = "\\HEEC, \(p\)-value" in 16; 
replace coeff = "HEEC, \(q\)-value" in 17; 

forvalues i = 1/9 { ;
	gen col`i' = "" ;
	label var col`i' "(`i')" ;
} ; 

* Fill in the columns ;
* Panel A ; 
local colnum = 1 ;
foreach y in score_tot score_cleanliness score_midwife score_condom { ;
	
	summ pval if outcome == "`y'" & coefficient == "_b[healthonly]" ;
	assert r(N) == 1 ;
	di r(mean) ;
	replace col`colnum' = string(r(mean), "%9.3f") in 2 ; 

	summ bh95_qval if outcome == "`y'" & coefficient == "_b[healthonly]" ;
	assert r(N) == 1 ;
	di r(mean) ;
	replace col`colnum' = string(r(mean), "%9.3f") in 3 ; 

	summ pval if outcome == "`y'" & coefficient == "_b[healthandpay]" ;
	assert r(N) == 1 ;
	replace col`colnum' = string(r(mean), "%9.3f") in 4 ; 
	
	summ bh95_qval if outcome == "`y'" & coefficient == "_b[healthandpay]" ;
	assert r(N) == 1 ;
	replace col`colnum' = string(r(mean), "%9.3f") in 5 ; 

	local ++colnum ;
} ;

* Panel B ; 
local colnum = 1 ;
foreach y in module_all module1_know module23_know module4_know module5_know { ;
	
	summ pval if outcome == "`y'" & coefficient == "_b[healthonly]" ;
	assert r(N) == 1 ;
	di r(mean) ;
	replace col`colnum' = string(r(mean), "%9.3f") in 9 ; 

	summ bh95_qval if outcome == "`y'" & coefficient == "_b[healthonly]" ;
	assert r(N) == 1 ;
	di r(mean) ;
	replace col`colnum' = string(r(mean), "%9.3f") in 10 ; 

	summ pval if outcome == "`y'" & coefficient == "_b[healthandpay]" ;
	assert r(N) == 1 ;
	replace col`colnum' = string(r(mean), "%9.3f") in 11 ; 
	
	summ bh95_qval if outcome == "`y'" & coefficient == "_b[healthandpay]" ;
	assert r(N) == 1 ;
	replace col`colnum' = string(r(mean), "%9.3f") in 12 ; 

	local ++colnum ;
} ;

* Panel C ; 
local colnum = 1 ;
foreach y in behav_all $module1_behav $module23_behav $module4_behav $module5_behav { ;
	
	summ pval if outcome == "`y'" & coefficient == "_b[healthonly]" ;
	assert r(N) == 1 ;
	di r(mean) ;
	replace col`colnum' = string(r(mean), "%9.3f") in 14 ; 

	summ bh95_qval if outcome == "`y'" & coefficient == "_b[healthonly]" ;
	assert r(N) == 1 ;
	di r(mean) ;
	replace col`colnum' = string(r(mean), "%9.3f") in 15 ; 

	summ pval if outcome == "`y'" & coefficient == "_b[healthandpay]" ;
	assert r(N) == 1 ;
	replace col`colnum' = string(r(mean), "%9.3f") in 16 ; 
	
	summ bh95_qval if outcome == "`y'" & coefficient == "_b[healthandpay]" ;
	assert r(N) == 1 ;
	replace col`colnum' = string(r(mean), "%9.3f") in 17 ; 

	local ++colnum ;
} ;

* output the table ;
keep coeff col1-col9 ;
dropmiss, obs force ;

texsave using "./output/appendix-table-multiple-inference-adjustment.tex",
	replace frag nofix
	title("Multiple Inference Adjustment")
	marker(table-multiple-inference-adjustment)
	varlabels
	align(lCCCCCCCCCC)
	width(0.85\linewidth)
	size(small)
	footnote("\textit{Notes:} 
	The table presents the adjusted and non-adjusted \(p\)-values of the coefficients 
	for HEE and HEEC.
	The adjusted \(p\)-values are the \textcite{benjamini1995controlling} \(q\)-values 
	as described in \textcite{anderson2008multiple}. Panels A, B, and C correspond to 
	Tables \ref{table-midline-knowledge-aggregate}, \ref{table-endline-knowledge-aggregate}, and
	\ref{table-endline-behavior}, respectively. For example, in Panel A, Column 1, the value 0.000 
	is the \(p\)-value of the estimate of the coefficient on HEE in 
	Table \ref{table-midline-knowledge-aggregate}, Column 1, and its corresponding \(q\)-value is 0.001.", size(footnotesize) width(p{0.85\linewidth})) ;

* cleanup table footer ;
filefilter "./output/appendix-table-multiple-inference-adjustment.tex" "./output/appendix-table-multiple-inference-adjustment-v00.tex", 
	from("\BSmulticolumn{10}{p{0.85\BSlinewidth}}{\BSbegin{footnotesize}") to("\BScaption*{\BSfootnotesize") replace ;
filefilter "./output/appendix-table-multiple-inference-adjustment-v00.tex" "./output/appendix-table-multiple-inference-adjustment-v01.tex", from("\BSend{footnotesize}}") to("}") replace ;
filefilter "./output/appendix-table-multiple-inference-adjustment-v01.tex" "./output/appendix-table-multiple-inference-adjustment-v02.tex", from("\BSend{tabularx}") to("") replace ;
filefilter "./output/appendix-table-multiple-inference-adjustment-v02.tex" "./output/appendix-table-multiple-inference-adjustment.tex", from("\BSbottomrule \BSaddlinespace[\BSbelowrulesep]") to("\BSbottomrule \BSend{tabularx} \BScaptionsetup{width=0.85\BSlinewidth}") replace ;

erase "temp.dta" ;
forvalues i = 0/2 { ;
	erase "./output/appendix-table-multiple-inference-adjustment-v0`i'.tex" ;
} ;

exit ; 
